---
title: "Interactions between tsetse endosymbiont and Glossina pallidipes Salivary Gland
  Hypertrophy Virus (GpSGHV) on heterologous host of Glossina Species"
author: "Guler Demirbas-Uzel1, Antonios A. Augustinos, Vangelis Doudoumis, Andrew G.
  Parker, Kostas Bourtzis, and Adly M. M. Abd-Alla1*"
date: "01/105/2021"
output: word_document
---

```{r}

setwd("C:/Users/abdallaa/OneDrive - IAEA/My/Guler_2019/virus_symbio_Nov_2020/Revision/raw_data_final")


library(ggplot2)
library(lattice)
library(gcookbook)
library(datasets)
library(MASS)
library(survival)
library(rmarkdown)
library(knitr)
library(coxme)
library(lme4)
library(nlme)
library(MuMIn)
#install.packages("ggthemes") # Install 
library(ggthemes) # Load

```

##Normality for fig1a, b and C

```{r}
#Normality for fig1a

fig1a <- read.csv("Fig_1a.csv",sep=",", row.names=NULL)
fig1a=na.omit(fig1a)
levels(fig1a$Time)
head(fig1a)

#R test normality for ALL Wiggleworthia data 
Box=boxcox(relative_density~Sex*Time, data=fig1a, lambda = seq(-2,2,0.1))
Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]
lambda = Cox2[1, "Box.x"] 
lambda

Wigglesworthia <- with(fig1a,relative_density)
trans <- with(fig1a,transformed)
library(rcompanion)

plotNormalHistogram(Wigglesworthia)
qqnorm(Wigglesworthia,ylab="relative_density")
qqline(Wigglesworthia, 
       col="red")

plotNormalHistogram(trans)
qqnorm(trans,ylab="relative_density")
qqline(trans, 
       col="red")

#select the model

model1 <- lm(transformed ~ Time + Species, data = fig1a)
model2 <- lm(transformed ~ Time + Species + Sex, data = fig1a)
model3 <- lm(transformed ~ Time*Species, data = fig1a)
model4 <- lm(transformed ~ Time*Species*Sex, data = fig1a)

AICc(model1,model2,model3, model4)
anova(model3, model4)
summary(model4)
#====================================================================================

#Normality for fig1b

fig1b <- read.csv("Fig_1b.csv",sep=",", row.names=NULL)
fig1b=na.omit(fig1b)
levels(fig1b$Time)
head(fig1b)

#R test normality for ALL Sodalis data 
Box=boxcox(relative_density~Sex*Time, data=fig1b, lambda = seq(-2,2,0.1))
Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]
lambda = Cox2[1, "Box.x"] 
lambda

Sodalis <- with(fig1b,relative_density)
trans <- with(fig1b,transformed)
library(rcompanion)

plotNormalHistogram(Sodalis)
qqnorm(Sodalis,ylab="relative_density")
qqline(Sodalis, 
       col="red")

plotNormalHistogram(trans)
qqnorm(trans,ylab="relative_density")
qqline(trans, 
       col="red")

#select the model

model1 <- lm(transformed ~ Time + Species, data = fig1b)
model2 <- lm(transformed ~ Time + Species + Sex, data = fig1b)
model3 <- lm(transformed ~ Time*Species, data = fig1b)
model4 <- lm(transformed ~ Time*Species*Sex, data = fig1b)

AICc(model1,model2,model3, model4)
anova(model3, model4)
summary(model4)

#==============================================================

#Normality for fig1c

fig1c <- read.csv("Fig_1c.csv",sep=",", row.names=NULL)
fig1c=na.omit(fig1c)
levels(fig1c$Time)
head(fig1c)

#R test normality for ALL Sodalis data 
Box=boxcox(relative_density~Sex*Time, data=fig1c, lambda = seq(-2,2,0.1))
Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]
lambda = Cox2[1, "Box.x"] 
lambda

Wolbachia <- with(fig1c,relative_density)
trans <- with(fig1c,transformed)
library(rcompanion)

plotNormalHistogram(Wolbachia)
qqnorm(Wolbachia,ylab="relative_density")
qqline(Wolbachia, 
       col="red")

plotNormalHistogram(trans)
qqnorm(trans,ylab="relative_density")
qqline(trans, 
       col="red")

#select the model

model1 <- lm(transformed ~ Time + Species, data = fig1c)
model2 <- lm(transformed ~ Time + Species + Sex, data = fig1c)
model3 <- lm(transformed ~ Time*Species, data = fig1c)
model4 <- lm(transformed ~ Time*Species*Sex, data = fig1c)

AICc(model1,model2,model3, model4)
anova(model3, model4)
summary(model4)
```

##Prepare Figure 1a, b and C

```{r}
fig1a <- read.csv("Fig_1a.csv",sep=",", row.names=NULL)
fig_1a<- ggplot(fig1a, aes(x=Time, y=transformed)) + 
  geom_point() + 
  geom_smooth(se = FALSE, method = lm) +  
  geom_vline(xintercept = 0.0, color = "black", size=1)+
  facet_grid(~Species) + ylim(-4,6.6)+ xlim(0,10)


tiff("fig_1a.tiff", width = 8, height = 4, units = 'in', res = 300)
plot(fig_1a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bold("Time post Injection (Days)"))) + ylab(expression(bolditalic("Wigglesworthia ")~bold("transformed relative density")))
dev.off()

fig1b <- read.csv("Fig_1b.csv",sep=",", row.names=NULL)
fig_1b<- ggplot(fig1b, aes(x=Time, y=transformed)) + 
  geom_point() + 
  geom_smooth(se = FALSE, method = lm) +  
  geom_vline(xintercept = 0.0, color = "black", size=1)+
  facet_grid(~Species) + ylim(-4,6.6)+ xlim(0,10)

tiff("fig_1b.tiff", width = 8, height = 4, units = 'in', res = 300)
plot(fig_1b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bold("Days post Injection (dpi)"))) + ylab(expression(bolditalic("Sodalis ")~bold("transformed relative density")))
dev.off()

fig1c <- read.csv("Fig_1c.csv",sep=",", row.names=NULL)
fig_1c<- ggplot(fig1c, aes(x=Time, y=transformed)) + 
  geom_point() + 
  geom_smooth(se = FALSE, method = lm) +  
  geom_vline(xintercept = 0.0, color = "black", size=1)+
  facet_grid(~Species) + ylim(-4,6.6)+ xlim(0,10)

tiff("fig_1c.tiff", width = 8, height = 4, units = 'in', res = 300)
plot(fig_1c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bold("Days post Injection (dpi)"))) + ylab(expression(bolditalic("Wolbachia ")~bold("transformed relative density")))
dev.off()

```

## #plotting supplementary figure 2a, b and C

```{r cars}
# Sup. Figure 1b, Line Wigglesworthia
(p <- ggplot(data = fig1a, aes(x = Sex, y = Species)) 
 + geom_tile(aes(fill = transformed), colour = "white")
 + scale_fill_gradient(low= "yellow", high ="black"))

# Sup. Figure 1a, SODALIS
(p <- ggplot(data = fig1b, aes(x =Sex, y = Species)) 
  + geom_tile(aes(fill = transformed), colour = "white")
  + scale_fill_gradient(low= "yellow", high ="black"))

# Sup. Figure 1c, WOlbachia 
(p <- ggplot(data = fig1c, aes(x = Sex, y = Species)) 
  + geom_tile(aes(fill = transformed), colour = "white")
  + scale_fill_gradient(low= "yellow", high ="black"))

```

# #Normality for supplementary figure 1 (figs1)

```{r}
figs1 <- read.csv("Fig_S1.csv",sep=",", row.names=NULL)
fig1c=na.omit(figs1)
levels(figs1$Time)
head(figs1)

#R test normality for ALL GpSGHV data 
Box=boxcox(relative_density~Sex*Time, data=fig1c, lambda = seq(-2,2,0.1))
Cox = data.frame(Box$x, Box$y)
Cox2 = Cox[with(Cox, order(-Cox$Box.y)),]
Cox2[1,]
lambda = Cox2[1, "Box.x"] 
lambda

GpSGHV <- with(figs1,relative_density)
trans <- with(figs1,transformed)
library(rcompanion)

plotNormalHistogram(GpSGHV)
qqnorm(GpSGHV,ylab="relative_density")
qqline(Wolbachia, 
       col="red")

plotNormalHistogram(trans)
qqnorm(trans,ylab="relative_density")
qqline(trans, 
       col="red")

#select the model

model1 <- lm(transformed ~ Time + Species, data = figs1)
model2 <- lm(transformed ~ Time + Species + Sex, data = figs1)
model3 <- lm(transformed ~ Time*Species, data = figs1)
model4 <- lm(transformed ~ Time*Species*Sex, data = figs1)

AICc(model1,model2,model3, model4)
anova(model3, model4)
summary(model4)

```

# #plotting supplementary 1 (figs1)

```{r}
figs1 <- read.csv("Fig_S1.csv",sep=",", row.names=NULL)
fig_s1<- ggplot(figs1, aes(x=Time, y=transformed)) + 
  geom_point() + 
  geom_smooth(se = FALSE, method = lm) +  
  geom_vline(xintercept = 0.0, color = "black", size=1)+
  facet_grid(~Species) + ylim(-4,6.6)+ xlim(0,10)
figs1

tiff("fig_s1.tiff", width = 8, height = 4, units = 'in', res = 300)
plot(fig_s1+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bold("Days post Injection (dpi)"))) + ylab(expression(bold("GpSGHV transformed relative density")))
dev.off()
```

## statistics for supplementary table 2

```{r}
###Statistics_Wigglesworthia_ANova test_Supplementary Table 2 

lm(formula = transformed ~ factor(Time)*Species, data = fig1a)

anova(with(fig1a,lm((transformed~Sex*Time*Species))))
anova(with(fig1a,lm((transformed~Time*Species))))  
anova(with(fig1a,lm(transformed~Sex*Species)))

#HSD analysis in case of significant differenes
fig1ahsd <-aov(with(fig1a,lm(transformed~factor(Time)*Species)))
TukeyHSD(fig1ahsd,ordered = FALSE,conf.level = 0.95)

lm(formula = transformed ~ Time, data = fig1a)
fit <-lm(transformed ~ Time, data = fig1a)
summary (fit)


###Statistics_Sodalis_ANova test_Supplementary Table 2 

lm(formula = transformed ~ factor(Time)*Species, data = fig1b)

anova(with(fig1b,lm((transformed~Sex*Time*Species))))
anova(with(fig1b,lm((transformed~Time*Species))))  
anova(with(fig1b,lm(transformed~Sex*Species)))

#HSD analysis in case of significant differenes
fig1bhsd <-aov(with(fig1b,lm(transformed~factor(Time)*Species)))
TukeyHSD(fig1bhsd,ordered = FALSE,conf.level = 0.95)

lm(formula = transformed ~ Time, data = fig1b)
fit <-lm(transformed ~ Time, data = fig1b)
summary (fit)

###Statistics_WOlbachia _ANova test_Supplementary Table 2 

lm(formula = transformed ~ factor(Time)*Species, data = fig1c)

anova(with(fig1c,lm((transformed~Sex*Time*Species))))
anova(with(fig1c,lm((transformed~Time*Species))))  
anova(with(fig1c,lm(transformed~Sex*Species)))

#HSD analysis in case of significant differences
fig1chsd <-aov(with(fig1c,lm(transformed~factor(Time)*Species)))
TukeyHSD(fig1chsd,ordered = FALSE,conf.level = 0.95)


lm(formula = transformed ~ Time, data = fig1c)
fit <-lm(transformed ~ Time, data = fig1c)
summary (fit)

###Statistics_GpSGHV denisty_Supplementary Figure 1 _ANova test_Supplementary Table 2 
fig1 <- read.csv("Fig_S1.csv",sep=",", row.names=NULL)
lm(formula = transformed ~ factor(Time)*Species, data = fig1)

anova(with(fig1,lm((transformed~Sex*Time*Species))))
anova(with(fig1,lm((transformed~Time*Species))))  
anova(with(fig1,lm(transformed~Sex*Species)))

#HSD analysis in case of significant differenes
fig1hsd <-aov(with(fig1,lm(transformed~factor(Time)*Species)))
TukeyHSD(fig1hsd,ordered = FALSE,conf.level = 0.95)


lm(formula = transformed ~ Time, data = fig1)
fit <-lm(transformed ~ Time, data = fig1)
summary (fit)

```

# #Statistics for supplementary table 3
```{r}
###Statistics_Wigglesworthia_ANova test(Female and Male)_Supplementary Table 3 

#Gb_Wigglesworthia_Anova_Female-Male
WigGb <- read.csv("Fig_1a.csv")
WigGb <- subset(WigGb, Species=="Gb")
anova(with(WigGb,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGb)


#Gf_Wigglesworthia_Anova_Female-Male
WigGff <- read.csv("Fig_1a.csv")
WigGff <- subset(WigGff, Species=="Gff")
anova(with(WigGff,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGff)


#Gmc_Wigglesworthia_Anova_Female-Male
WigGmc<- read.csv("Fig_1a.csv")
WigGmc <- subset(WigGmc, Species=="Gmc")
anova(with(WigGmc,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGmc)

#Gmm_Wigglesworthia_Anova_Female-Male
WigGmm<- read.csv("Fig_1a.csv")
WigGmm <- subset(WigGmm, Species=="Gmm")
anova(with(WigGmm,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGmm)


#Gp_Wigglesworthia_Anova_Female-Male
WigGp<- read.csv("Fig_1a.csv")
WigGp <- subset(WigGp, Species=="Gp")
anova(with(WigGp,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGp)

#Gpg_Wigglesworthia_Anova_Female-Male
WigGpg<- read.csv("Fig_1a.csv")
WigGpg <- subset(WigGpg, Species=="Gpg")
anova(with(WigGpg,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WigGpg)

###Statistics_Sodalis_ANova test(Female and Male)_Supplementary Table 3 

#Gb_Sodalis_Anova_Female-Male
SodGb <- read.csv("Fig_1b.csv")
SodGb <- subset(SodGb, Species=="Gb")
anova(with(SodGb,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGb)

#Gf_Sodalis_Anova_Female-Male
SodGff <- read.csv("Fig_1b.csv")
SodGff <- subset(SodGff, Species=="Gf")
anova(with(SodGff,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGff)

#Gmc_Sodalis_Anova_Female-Male
SodGmc<- read.csv("Fig_1b.csv")
SodGmc <- subset(SodGmc, Species=="Gmc")
anova(with(SodGmc,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGmc)

#Gmm_Sodalis_Anova_Female-Male
SodGmm<- read.csv("Fig_1b.csv")
SodGmm <- subset(SodGmm, Species=="Gmm")
anova(with(WigGmm,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGmm)

#Gp_Sodalis_Anova_Female-Male
SodGp<- read.csv("Fig_1b.csv")
SodGp <- subset(SodGp, Species=="Gp")
anova(with(SodGp,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGp)

#Gpg_Sodalis_Anova_Female-Male
SodGpg<- read.csv("Fig_1b.csv")
SodGpg <- subset(SodGpg, Species=="Gpg")
anova(with(SodGpg,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = SodGpg)

###Statistics_Wolbachia_ANova test(Female and Male)_Supplementary Table 3 

#Gb_Wolbachia_Anova_Female-Male
WolGb <- read.csv("Fig_1c.csv")
WolGb <- subset(WolGb, Species=="Gb")
anova(with(WolGb,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGb)

#Gf_Wolbachia_Anova_Female-Male
WolGff <- read.csv("Fig_1c.csv")
WolGff <- subset(WolGff, Species=="Gff")
anova(with(WolGff,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGff)

#Gmc_Wolbachia_Anova_Female-Male
WolGmc<- read.csv("Fig_1c.csv")
WolGmc <- subset(WolGmc, Species=="Gmc")
anova(with(WolGmc,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGmc)

#Gmm_Wolbachia_Anova_Female-Male
WolGmm<- read.csv("Fig_1c.csv")
WolGmm <- subset(WolGmm, Species=="Gmm")
anova(with(WolGmm,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGmm)

#Gp_Wolbachia_Anova_Female-Male
WolGp<- read.csv("Fig_1c.csv")
WolGp <- subset(WolGp, Species=="Gp")
anova(with(WolGp,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGp)

#Gpg_Wolbachia_Anova_Female-Male
WolGpg<- read.csv("Fig_1c.csv")
WolGpg <- subset(WolGpg, Species=="Gpg")
anova(with(SodGp,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = WolGpg)

###Statistics_GpSGHV_Anova test_(Female and Male)_Supplementary Table 3

#Gb_GpSGHV_Anova_Female-Male
GpSGHVGb <- read.csv("Fig_S1.csv")
GpSGHVGb <- subset(GpSGHVGb, Species=="Gb")
anova(with(GpSGHVGb,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGb)

#Gf_GpSGHV_Anova_Female-Male
GpSGHVGff <- read.csv("Fig_S1.csv")
GpSGHVGff <- subset(GpSGHVGff, Species=="Gff")
anova(with(GpSGHVGff,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGff)

#Gmc_GpSGHV_Anova_Female-Male
GpSGHVGmc<- read.csv("Fig_S1.csv")
GpSGHVGmc <- subset(GpSGHVGmc, Species=="Gmc")
anova(with(GpSGHVGmc,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGmc)

#Gmm_GpSGHV_Anova_Female-Male
GpSGHVGmm<- read.csv("Fig_S1.csv")
GpSGHVGmm <- subset(GpSGHVGmm, Species=="Gmm")
anova(with(GpSGHVGmm,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGmm)

#Gp_GpSGHV_Anova_Female-Male
GpSGHVGp<- read.csv("Fig_S1.csv")
GpSGHVGp <- subset(GpSGHVGp, Species=="Gp")
anova(with(GpSGHVGp,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGp)

#Gpg_GpSGHV_Anova_Female-Male
GpSGHVGpg<- read.csv("Fig_S1.csv")
GpSGHVGpg <- subset(GpSGHVGpg, Species=="Gpg")
anova(with(GpSGHVGpg,lm(transformed~Sex*Time)))
lm(formula = transformed ~ Time, data = GpSGHVGpg)


```

# Statistics for supplementary table 4

```{r}
#Statistics for supplementary table 4
#Wigglesworthia Regression

WigGmm <- read.csv("fig_1a.csv")
WigGmm <- subset(WigGmm, Species=="Gmm")
lm(formula = transformed ~ Time, data = WigGmm)
fit <-lm(transformed ~ Time, data = WigGmm)
summary (fit)

###Wigglesworthia Gmc_Regression

WigGmc <- read.csv("fig_1a.csv")
WigGmc <- subset(WigGmc, Species=="Gmc")
lm(formula = transformed ~ Time, data = WigGmc)
fit <-lm(transformed ~ Time, data = WigGmc)
summary (fit)

###Wigglesworthia Gb_Regression

WigGb <- read.csv("fig_1a.csv")
WigGb <- subset(WigGb, Species=="Gb")
lm(formula = transformed ~ Time, data = WigGb)
fit <-lm(transformed ~ Time, data = WigGb)
summary (fit)
###Wigglesworthia Gpg_Regression

WigGpg <- read.csv("fig_1a.csv")
WigGpg <- subset(WigGpg, Species=="Gpg")
lm(formula = transformed ~ Time, data = WigGpg)
fit <-lm(transformed ~ Time, data = WigGpg)
summary (fit)
###Wigglesworthia Gff_Regression

WigGff <- read.csv("fig_1a.csv")
WigGff <- subset(WigGff, Species=="Gff")
lm(formula = transformed ~ Time, data = WigGff)
fit <-lm(transformed ~ Time, data = WigGff)
summary (fit)
###Wigglesworthia Gp_Regression

WigGp <- read.csv("fig_1a.csv")
WigGp <- subset(WigGp, Species=="Gp")
lm(formula = transformed ~ Time, data = WigGp)
fit <-lm(transformed ~ Time, data = WigGp)
summary (fit)
#-----------------------------------------------
#Sodalis Regression

###Sodalis Gmm_Regression

SodGmm <- read.csv("Fig_1b.csv")
SodGmm <- subset(SodGmm, Species=="Gmm")
lm(formula = transformed ~ Time, data = SodGmm)
fit <-lm(transformed ~ Time, data = SodGmm)
summary (fit)

###Sodalis Gmc_Regression

SodGmc <- read.csv("Fig_1b.csv")
SodGmc <- subset(SodGmc, Species=="Gmc")
lm(formula = transformed ~ Time, data = SodGmc)
fit <-lm(transformed ~ Time, data = SodGmc)
summary (fit)

###Sodalis Gb_Regression

SodGb <- read.csv("Fig_1b.csv")
SodGb <- subset(SodGb, Species=="Gb")
lm(formula = transformed ~ Time, data = SodGb)
fit <-lm(transformed ~ Time, data = SodGb)
summary (fit)
###Sodalis Gpg_Regression

SodGpg <- read.csv("Fig_1b.csv")
SodGpg <- subset(SodGpg, Species=="Gpg")
lm(formula = transformed ~ Time, data = SodGpg)
fit <-lm(transformed ~ Time, data = SodGpg)
summary (fit)
###Sodalis Gff_Regression

SodGf <- read.csv("Fig_1b.csv")
SodGf <- subset(SodGf, Species=="Gf")
lm(formula = transformed ~ Time, data = SodGf)
fit <-lm(transformed ~ Time, data = SodGf)
summary (fit)
###Sodalis Gp_Regression

SodGp <- read.csv("Fig_1b.csv")
SodGp <- subset(SodGp, Species=="Gp")
lm(formula = transformed ~ Time, data = SodGp)
fit <-lm(transformed ~ Time, data = SodGp)
summary (fit)

#----------------------------------------------------------
#Wolbachia regression


###WOlbachia  Gmm_Regression

WolGmm <- read.csv("Fig_1c.csv")
WolGmm <- subset(WolGmm, Species=="Gmm")
lm(formula = transformed ~ Time, data = WolGmm)
fit <-lm(transformed ~ Time, data = WolGmm)
summary (fit)

### WOlbachia  Gmc_Regression

WolGmc <- read.csv("Fig_1c.csv")
WolGmc <- subset(WolGmc, Species=="Gmc")
lm(formula = transformed ~ Time, data = WolGmc)
fit <-lm(transformed ~ Time, data = WolGmc)
summary (fit)

### WOlbachia  Gb_Regression

WolGb <- read.csv("Fig_1c.csv")
WolGb <- subset(WolGb, Species=="Gb")
lm(formula = transformed ~ Time, data = WolGb)
fit <-lm(transformed ~ Time, data = WolGb)
summary (fit)
### WOlbachia  Gpg_Regression

WolGpg <- read.csv("Fig_1c.csv")
WolGpg <- subset(WolGpg, Species=="Gpg")
lm(formula = transformed ~ Time, data = WolGpg)
fit <-lm(transformed ~ Time, data = WolGpg)
summary (fit)
### WOlbachia  Gff_Regression

WolGff <- read.csv("Fig_1c.csv")
WolGff <- subset(WolGff, Species=="Gff")
lm(formula = transformed ~ Time, data = WolGff)
fit <-lm(transformed ~ Time, data = WolGff)
summary (fit)
### WOlbachia  Gp_Regression

WolGp <- read.csv("Fig_1c.csv")
WolGp <- subset(WolGp, Species=="Gp")
lm(formula = transformed ~ Time, data = WolGp)
fit <-lm(transformed ~ Time, data = WolGp)
summary (fit)

#-----------------------------------------------------------------

# GpSGHV Regression

wolGmc <- read.csv("Fig_S1.csv")
WolGmc <- subset(WolGmc, Species=="Gmc")
lm(formula = transformed ~ Time, data = WolGmc)
fit <-lm(transformed ~ Time, data = WolGmc)
summary (fit)

### GpSGHV  Gb_Regression

WolGb <- read.csv("Fig_S1.csv")
WolGb <- subset(WolGb, Species=="Gb")
lm(formula = transformed ~ Time, data = WolGb)
fit <-lm(transformed ~ Time, data = WolGb)
summary (fit)


### GpSGHV  Gpg_Regression

WolGpg <- read.csv("Fig_S1.csv")
WolGpg <- subset(WolGpg, Species=="Gpg")
lm(formula = transformed ~ Time, data = WolGpg)
fit <-lm(transformed ~ Time, data = WolGpg)
summary (fit)

### GpSGHV  Gff_Regression

WolGff <- read.csv("Fig_S1.csv")
WolGff <- subset(WolGff, Species=="Gff")
lm(formula = transformed ~ Time, data = WolGff)
fit <-lm(transformed ~ Time, data = WolGff)
summary (fit)

### WOlbachia  Gp_Regression

WolGp <- read.csv("Fig_S1.csv")
WolGp <- subset(WolGp, Species=="Gp")
lm(formula = transformed ~ Time, data = WolGp)
fit <-lm(transformed ~ Time, data = WolGp)
summary (fit)

```

# Plotting Figure 2

```{r}
#plotting fig 2

fig2 <- read.csv("Fig_2.csv")
fig2a<-ggplot(fig2, aes(Wigglesworthia, GpSGHV, col = Species, fill = Species)) + 
  geom_point(size = 3, shape = 21, col = "black") +
  geom_vline(xintercept = 2.1, color = "black", size=1)+
  xlab(expression(italic("Wigglesworthis")))

tiff("fig2a.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig2a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Wigglesworthia"))) + ylab(expression(bold("GpSGHV")))
dev.off()

fig2b<-ggplot(fig2, aes(Sodalis, GpSGHV, col = Species, fill = Species)) + 
  geom_point(size = 3, shape = 21, col = "black") +
  geom_vline(xintercept = 18, color = "black", size=1)

tiff("fig2b.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig2b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Sodalis"))) + ylab(expression(bold("GpSGHV")))
dev.off()

fig2c<-ggplot(fig2, aes(Wolbachia, GpSGHV, col = Species, fill = Species)) + 
  geom_point(size = 3, shape = 21, col = "black") +
  geom_vline(xintercept = 18, color = "black", size=1)

tiff("fig2c.tiff", width = 4, height = 4, units = 'in', res = 300)
plot(fig2c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black"))+ theme(legend.position = c(.95, .95),legend.justification = c("right", "top"))) + xlab(expression(bolditalic("Sodalis"))) + ylab(expression(bold("GpSGHV")))
dev.off()


```

